#############################################################################################
# summarize NKPH estimates
#############################################################################################

# activate project environment
cd( dirname(Base.source_path()) )
using Pkg
Pkg.activate("../")

using ArgParse
using DataFrames
using XLSX

println("loading NKPH methods")
include("../src/NewKeynesianPreferredHabitatModels.jl");
using .NewKeynesianPreferredHabitatModels


# parse command line arguments
s = ArgParseSettings()
@add_arg_table! s begin
    # input options
    "--run_summary"
        help = "summarize model"
        arg_type = Int
        default = 1
    "--run_baseline"
        help = "baseline QE/QT experiments"
        arg_type = Int
        default = 1
    "--run_robustness"
        help = "robustness QE/QT experiments"
        arg_type = Int
        default = 1
    "--N_param_vals_robustness"
        help = "number of parameter runs for each robustness experiment"
        arg_type = Int
        default = 31
end




#############################################################################################
# helper functions for robustness exercises
function _reset_nkph_param!(
        param_name::String,
        param_type::Symbol,
        param_val::Float64,
        nkphopt::NKPHEstimationOptimizer
    )
    # update specified param
    if param_type == :fix
        idx_param = findfirst(nkphopt.df_params_fixed[:, :param] .== param_name)
        nkphopt.df_params_fixed[idx_param, :val] = param_val
    elseif param_type == :est
        idx_param = findfirst(nkphopt.df_params_constraints[:, :param] .== param_name)
        nkphopt.df_params_constraints[idx_param, :init] = param_val
    elseif param_type == :risk_aversion
        idx_param = findfirst(nkphopt.df_soln_opts[:, :continuation_opts] .== "tracking_target")
        nkphopt.df_soln_opts[idx_param, :value] = param_val
    else
        throw(DomainError(param_type, "param_type must be fix/est/risk_aversion"))
    end
    return nothing
end

function reload_nkph_model_estimates!(
        param_name::String,
        param_type::Symbol,
        param_val::Float64,
        nkphopt::NKPHEstimationOptimizer
    )
    # update specified (single) param and re-solve
    _reset_nkph_param!(param_name, param_type, param_val, nkphopt)
    _init_nkphopt(nkphopt)
    return nothing
end

function reload_nkph_model_estimates!(
        param_names::Vector{String},
        param_types::Vector{Symbol},
        param_vals::Vector{Float64},
        nkphopt::NKPHEstimationOptimizer
    )
    # update specified (multiple) params and re-solve
    for (param_name, param_type, param_val) in zip(param_names, param_types, param_vals)
        _reset_nkph_param!(param_name, param_type, param_val, nkphopt)
    end
    _init_nkphopt(nkphopt)
    return nothing
end

#############################################################################################
# size of auction/QE shock function

function estimate_b_shock_size(
        short_idx::Int, 
        long_idx::Int, 
        sig_bs_scale::Float64,
        nkphopt::NKPHEstimationOptimizer
    )::Tuple{Float64, Float64}
    # estimate implied shock sizes from alt localization coeffs
    # pull out objects
    tau_arr, tau_arr_reg = nkphopt.tau_arr, nkphopt.tau_arr_reg
    A_tau_arr = nkphopt.nkphmmts.A_tau_arr
    localization_short_hat = nkphopt.nkphemp.df_alt_regs[:, :alt_localization_short]
    localization_long_hat = nkphopt.nkphemp.df_alt_regs[:, :alt_localization_long]

    # data and model-implied value (using correct tau index)
    numer_val = 0.
    denom_val = 0.
    loc_idx = 0
    for (i, tau) in enumerate(tau_arr)
        if tau in tau_arr_reg
            loc_idx += 1
            loc_short = localization_short_hat[loc_idx]
            loc_long = localization_long_hat[loc_idx]

            A_tau_scaled = 100 * A_tau_arr[:, i] / tau
            numer_val += sig_bs_scale * A_tau_scaled[short_idx] * loc_short
            numer_val += A_tau_scaled[long_idx] * loc_long

            denom_val += (sig_bs_scale * A_tau_scaled[short_idx])^2
            denom_val += (A_tau_scaled[long_idx])^2
        end
    end
    sig_bl = numer_val / denom_val
    sig_bs = sig_bs_scale * sig_bl
    return sig_bl, sig_bs
end

function estimate_sig_auc_size(nkphopt::NKPHEstimationOptimizer;
        sig_bs_scale=1.5, short_idx=5, long_idx=6
    )
    # compute implied size given auction shock info
    return estimate_b_shock_size(short_idx, long_idx, sig_bs_scale, nkphopt)
end

function estimate_sig_QE_size(nkphopt::NKPHEstimationOptimizer;
        sig_bs_scale=1.5, short_idx=5, long_idx=6,
        sig_QE_scale=50.
    )
    # compute implied size of QE given auction shock info
    sig_bl, sig_bs = estimate_sig_auc_size(nkphopt;
        sig_bs_scale=sig_bs_scale, short_idx=short_idx, long_idx=long_idx
    )
    sig_QE = -sig_QE_scale * sig_bl
    #println(sig_bl, " ", sig_bs, " ", sig_QE)
    return sig_QE
end


#############################################################################################
# QE experiments

function run_experiment_baseline(
        shock_arr::Vector{Float64},
        nkphopt::NKPHEstimationOptimizer;
        t_arr=nothing, tau_arr=nothing, tau_val=10.
    )
    # macro IRF and yield curve responses to shock
    if t_arr === nothing
        t_arr = Vector{Float64}(range(0, 10., step=0.1))
    end
    if tau_arr === nothing
        tau_arr = Vector{Float64}(range(0.1, 30., step=0.1))
    end

    # yield responses (full yield curve, on impact)
    yc_resp, yc_resp_til = calc_yield_response(shock_arr, tau_arr, nkphopt.nkphsolvr)
    yc_mat = hcat(tau_arr, yc_resp, yc_resp_til)
    df_yc = DataFrame(yc_mat, [:tau, :y_tau, :ytil_tau])

    # yield responses (specified maturity, IRF)
    yc_resp_irf, yc_resp_til_irf = calc_yield_response_irf(shock_arr,
        t_arr, tau_val, nkphopt.nkphsolvr)
    yc_irf_mat = hcat(t_arr, yc_resp_irf, yc_resp_til_irf)
    df_yc_irf = DataFrame(yc_irf_mat, [:t, :y_tau, :ytil_tau])

    # IRF
    yt_irf, xt_irf = calc_irf(shock_arr, t_arr, nkphopt.nkphsolvr)
    yt_irf_mat = hcat(t_arr, yt_irf')
    xt_irf_mat = hcat(t_arr, xt_irf')
    # to dataframe
    yt_colnames = Symbol.(nkphopt.df_state_vars[:, :name][1:nkphopt.nkphsolvr.J])
    pushfirst!(yt_colnames, :t)
    xt_colnames = Symbol.(nkphopt.df_state_vars[:, :name][nkphopt.nkphsolvr.J+1:end])
    pushfirst!(xt_colnames, :t)
    df_yt_irf = DataFrame(yt_irf_mat, yt_colnames)
    df_xt_irf = DataFrame(xt_irf_mat, xt_colnames)
    return df_yc, df_yc_irf, df_yt_irf, df_xt_irf
end

##########################################
# robustness experiments
function run_experiment_robustness(
        param_names::Union{ String, Vector{String} },
        param_types::Union{ Symbol, Vector{Symbol} },
        param_vals_mat::Union{ Vector{Float64}, Matrix{Float64} },
        shock_arr::Vector{Float64},
        nkphopt::NKPHEstimationOptimizer;
        tau_resp=10., display=true
    )
    # robustness exercise, varying specified (multiple) params
    # save outcomes in matrix
    if display
        println("experiment: ", param_names)
    end
    N_param_vals = size(param_vals_mat, 1)
    if typeof(param_names)==String
        N_params = 1
    else
        N_params = length(param_names)
    end
    N_jumpvars = nkphopt.nkphsolvr.N - nkphopt.nkphsolvr.J
    outcome_mat = zeros(N_param_vals, N_params + 2 + 2*N_jumpvars)
    for idx=1:N_param_vals
        if display
            print(".")
        end
        if N_params==1
            param_vals = param_vals_mat[idx]
        else
            param_vals = param_vals_mat[idx, :]
        end
        # reload model and save outcomes
        reload_nkph_model_estimates!(param_names, param_types, param_vals, nkphopt)
        # yield response
        yc_resp, yc_resp_til = calc_yield_response(shock_arr, [tau_resp], nkphopt.nkphsolvr)
        # longrun response
        _, x_oo = calc_longrun_effect(shock_arr, nkphopt.nkphsolvr)
        # longrun volatility
        _, sd_x_oo = calc_longrun_volatilty(nkphopt.nkphmmts)
        # outcomes into vector
        outcome_arr = vcat( [yc_resp[1], yc_resp_til[1]], x_oo, sd_x_oo)
        # save to combined matrix
        if N_params==1
            outcome_mat[idx, 1] = param_vals
        else
            outcome_mat[idx, 1:N_params] = param_vals
        end
        outcome_mat[idx, N_params+1:end] = outcome_arr
    end
    if display
        println("")
    end
    # outcomes into dataframe
    xt_names = nkphopt.df_state_vars[:, :name][nkphopt.nkphsolvr.J+1:end]
    xt_oo_colnames = Symbol.( [xtn * "_oo" for xtn in xt_names] )
    sd_xt_oo_colnames = Symbol.( ["sd_" * xtn * "_oo" for xtn in xt_names] )
    colnames = vcat( Symbol.(param_names), [:y_tau, :ytil_tau],
        xt_oo_colnames, sd_xt_oo_colnames
    )
    return DataFrame(outcome_mat, colnames)
end


#############################################################################
# wrappers
function run_mp_experiments(mp_model_fname::String, xls_output_fname::String;
        mp_shock_size=-0.50, mp_shock_idx=nothing
    )
    # load MP model and run MP experiments
    # all experiments saved to xls file
    xls_output_fpath = OUTPUT_DIR * xls_output_fname * ".xlsx"
    # load model
    nkphopt_mp = load_nkph_model_estimates(mp_model_fname, true)
    if mp_shock_idx===nothing
        mp_shock_idx = nkphopt_mp.nkphsolvr.J
    end
    mp_shock = zeros(nkphopt_mp.nkphsolvr.J)
    mp_shock[mp_shock_idx] = mp_shock_size


    # create xls file, save basic info
    df_mp_summary = DataFrame(
        description = ["root_max", "mp_model_fname"],
        values = [maximum(abs.(nkphopt_mp.nkphsolvr.root_vec)), mp_model_fname]
    )
    # mp experiment
    df_mp_yc, df_mp_yc_irf, df_mp_yt_irf, df_mp_xt_irf =
        run_experiment_baseline(mp_shock, nkphopt_mp)
    # save to xls file
    XLSX.writetable(xls_output_fpath, overwrite=true, 
        "summary" => df_mp_summary, 
        "mp_yc" => df_mp_yc, 
        "mp_yc_irf" => df_mp_yc_irf, 
        "mp_yt_irf" => df_mp_yt_irf, 
        "mp_xt_irf" => df_mp_xt_irf, 
    )
    return nothing
end


function run_qe_experiments_all(qe_model_fname::String, xls_output_fname::String;
        run_baseline=true, run_robustness=false, run_macro=false,
        N_param_vals_robustness=51, tau_robustness=10.,
        qe_shock_size=nothing, qe_shock_idx=nothing,
        mp_shock_size=-0.25, mp_shock_idx=1
    )
    # load QE model and run all QE experiments
    # all experiments saved to xls file
    xls_output_fpath = OUTPUT_DIR * xls_output_fname * ".xlsx"
    # load model
    nkphopt_qe = load_nkph_model_estimates(qe_model_fname, true)
    # pull out default parameterization
    df_soln_opts_orig = copy(nkphopt_qe.df_soln_opts)
    df_params_constraints_orig = copy(nkphopt_qe.df_params_constraints)
    df_params_fixed_orig = copy(nkphopt_qe.df_params_fixed)

    function _reset_params_orig()
        # reset to original parameters before each experiment
        nkphopt_qe.df_soln_opts = copy(df_soln_opts_orig)
        nkphopt_qe.df_params_constraints = copy(df_params_constraints_orig)
        nkphopt_qe.df_params_fixed = copy(df_params_fixed_orig)
    end

    # construct qe shock
    qe_shock = zeros(nkphopt_qe.nkphsolvr.J)
    if qe_shock_size===nothing
        qe_shock_size = estimate_sig_QE_size(nkphopt_qe)
    end
    if qe_shock_idx===nothing
        qe_shock_idx = nkphopt_qe.nkphsolvr.J
    end
    qe_shock[qe_shock_idx] = qe_shock_size
    # comparison: mp shock
    mp_shock = zeros(nkphopt_qe.nkphsolvr.J)
    mp_shock[mp_shock_idx] = mp_shock_size


    if run_baseline
        # create xls file, save basic info
        df_qe_summary = DataFrame(
            description = ["root_max", "qe_model_fname"],
            values = [maximum(abs.(nkphopt_qe.nkphsolvr.root_vec)), qe_model_fname]
        )

        # baseline qe and mp experiment
        df_qe_yc, df_qe_yc_irf, df_qe_yt_irf, df_qe_xt_irf =
            run_experiment_baseline(qe_shock, nkphopt_qe)
        df_mp_yc, df_mp_yc_irf, df_mp_yt_irf, df_mp_xt_irf =
            run_experiment_baseline(mp_shock, nkphopt_qe)
        # save to xls file
        XLSX.writetable(xls_output_fpath, overwrite=true, 
            "summary" => df_qe_summary, 
            "qe_yc" => df_qe_yc, 
            "qe_yc_irf" => df_qe_yc_irf, 
            "qe_yt_irf" => df_qe_yt_irf, 
            "qe_xt_irf" => df_qe_xt_irf, 
            "mp_yc" => df_mp_yc, 
            "mp_yc_irf" => df_mp_yc_irf, 
            "mp_yt_irf" => df_mp_yt_irf, 
            "mp_xt_irf" => df_mp_xt_irf, 
        )
    end

    if run_robustness
        # mean reversion
        _reset_params_orig()
        experiment_name = "mean_reversion"
        param_name = "kappa_QE"
        param_type = :fix
        param_vals = Vector{Float64}(range(0.05, 1.0, length=N_param_vals_robustness))
        df_qe_experiment = run_experiment_robustness(
            param_name, param_type, param_vals, qe_shock, nkphopt_qe; tau_resp=tau_robustness
        )
        _df2xls_sheet!(df_qe_experiment, xls_output_fpath, "qe_exp_" * experiment_name)


        # QE volatility risk
        _reset_params_orig()
        experiment_name = "policy_volatility"
        param_name = "sigma_QE"
        param_type = :fix
        param_vals = Vector{Float64}(range(0.0, 50., length=N_param_vals_robustness))
        df_qe_experiment = run_experiment_robustness(
            param_name, param_type, param_vals, qe_shock, nkphopt_qe; tau_resp=tau_robustness
        )
        _df2xls_sheet!(df_qe_experiment, xls_output_fpath, "qe_exp_" * experiment_name)


        # QE maturity duration (note: multiple parameters)
        _reset_params_orig()
        experiment_name = "maturity_dur"
        param_names = ["theta1_QE", "theta1_til_QE"]
        param_types = [:fix, :fix]
        param_vals = hcat(
            Vector{Float64}(range(0.1, 2., length=N_param_vals_robustness)),
            Vector{Float64}(range(0.1, 2., length=N_param_vals_robustness))
        )
        df_qe_experiment = run_experiment_robustness(
            param_names, param_types, param_vals, qe_shock, nkphopt_qe; tau_resp=tau_robustness
        )
        _df2xls_sheet!(df_qe_experiment, xls_output_fpath, "qe_exp_" * experiment_name)


        # QE active/passive (note: multiple parameters)
        _reset_params_orig()
        experiment_name = "active_passive"
        param_names = ["theta0_QE", "theta0_til_QE", "theta0_QE_p", "theta0_til_QE_p"]
        param_types = [:fix, :fix, :fix, :fix, :fix]
        # scale weights from active to passive
        _scale_size_arr = Vector{Float64}(range(1., 0., length=N_param_vals_robustness))
        # default active magnitudes
        _theta0_QE_val = nkphopt_qe.df_params_fixed[
            findfirst(nkphopt_qe.df_params_fixed[:, :param] .== "theta0_QE"), :val]
        _theta0_til_QE_val = nkphopt_qe.df_params_fixed[
            findfirst(nkphopt_qe.df_params_fixed[:, :param] .== "theta0_til_QE"), :val]
        _param_theta0_QE_vals = _theta0_QE_val .* _scale_size_arr
        _param_theta0_til_QE_vals = _theta0_til_QE_val .* _scale_size_arr
        # scale down active magnitudes, and scale up passive magnitudes
        param_vals = hcat(
            _param_theta0_QE_vals, _param_theta0_til_QE_vals,
            reverse(_param_theta0_QE_vals), reverse(_param_theta0_til_QE_vals)
        )
        df_qe_experiment = run_experiment_robustness(
            param_names, param_types, param_vals, qe_shock, nkphopt_qe; tau_resp=tau_robustness
        )
        _df2xls_sheet!(df_qe_experiment, xls_output_fpath, "qe_exp_" * experiment_name)


        # risk aversion
        _reset_params_orig()
        experiment_name = "risk_aversion"
        param_name = "a"
        param_type = :risk_aversion
        param_vals = Vector{Float64}(range(0.025, 0.175, length=N_param_vals_robustness))
        df_qe_experiment = run_experiment_robustness(
            param_name, param_type, param_vals, qe_shock, nkphopt_qe; tau_resp=tau_robustness
        )
        _df2xls_sheet!(df_qe_experiment, xls_output_fpath, "qe_exp_" * experiment_name)
    end

    if run_macro
        # payoff risk
        _reset_params_orig()
        experiment_name = "payoff_risk"
        param_name = "s_d"
        param_type = :est
        param_vals = Vector{Float64}(range(0.5, 3., length=N_param_vals_robustness))
        df_qe_experiment = run_experiment_robustness(
            param_name, param_type, param_vals, qe_shock, nkphopt_qe; tau_resp=tau_robustness
        )
        _df2xls_sheet!(df_qe_experiment, xls_output_fpath, "qe_exp_" * experiment_name)


        # effective borrowing rate shape
        _reset_params_orig()
        experiment_name = "effrate_shape"
        param_name = "eta1_til"
        param_type = :fix
        param_vals = Vector{Float64}(range(0.5, 4., length=N_param_vals_robustness))
        df_qe_experiment = run_experiment_robustness(
            param_name, param_type, param_vals, qe_shock, nkphopt_qe; tau_resp=tau_robustness
        )
        _df2xls_sheet!(df_qe_experiment, xls_output_fpath, "qe_exp_" * experiment_name)


        # effective borrowing rate weight
        _reset_params_orig()
        experiment_name = "effrate_weight"
        param_name = "n0"
        param_type = :est
        param_vals = Vector{Float64}(range(0.0, 1.0, length=N_param_vals_robustness))
        df_qe_experiment = run_experiment_robustness(
            param_name, param_type, param_vals, qe_shock, nkphopt_qe; tau_resp=tau_robustness
        )
        _df2xls_sheet!(df_qe_experiment, xls_output_fpath, "qe_exp_" * experiment_name)


        # taylor rule coeff: pi
        _reset_params_orig()
        experiment_name = "mp_pi"
        param_name = "p_ipi"
        param_type = :est
        param_vals = Vector{Float64}(range(1.5, 4., length=N_param_vals_robustness))
        df_qe_experiment = run_experiment_robustness(
            param_name, param_type, param_vals, qe_shock, nkphopt_qe; tau_resp=tau_robustness
        )
        _df2xls_sheet!(df_qe_experiment, xls_output_fpath, "qe_exp_" * experiment_name)


        # taylor rule coeff: x
        _reset_params_orig()
        experiment_name = "mp_x"
        param_name = "p_ix"
        param_type = :est
        param_vals = Vector{Float64}(range(0.0, 1., length=N_param_vals_robustness))
        df_qe_experiment = run_experiment_robustness(
            param_name, param_type, param_vals, qe_shock, nkphopt_qe; tau_resp=tau_robustness
        )
        _df2xls_sheet!(df_qe_experiment, xls_output_fpath, "qe_exp_" * experiment_name)


        # taylor rule coeff: mean reversion
        _reset_params_orig()
        experiment_name = "mp_reversion"
        param_name = "k_i"
        param_type = :est
        param_vals = Vector{Float64}(range(0.4, 2.5, length=N_param_vals_robustness))
        df_qe_experiment = run_experiment_robustness(
            param_name, param_type, param_vals, qe_shock, nkphopt_qe; tau_resp=tau_robustness
        )
        _df2xls_sheet!(df_qe_experiment, xls_output_fpath, "qe_exp_" * experiment_name)


        # rigidity coeff
        _reset_params_orig()
        experiment_name = "pc_rigidity"
        param_name = "delta_x"
        param_type = :est
        param_vals = Vector{Float64}(range(0.4, 1.0, length=N_param_vals_robustness))
        df_qe_experiment = run_experiment_robustness(
            param_name, param_type, param_vals, qe_shock, nkphopt_qe; tau_resp=tau_robustness
        )
        _df2xls_sheet!(df_qe_experiment, xls_output_fpath, "qe_exp_" * experiment_name)


        # dividend process coeff: output
        _reset_params_orig()
        experiment_name = "dividend_x"
        param_name = "p_dx"
        param_type = :est
        param_vals = Vector{Float64}(range(0.0, 1., length=N_param_vals_robustness))
        df_qe_experiment = run_experiment_robustness(
            param_name, param_type, param_vals, qe_shock, nkphopt_qe; tau_resp=tau_robustness
        )
        _df2xls_sheet!(df_qe_experiment, xls_output_fpath, "qe_exp_" * experiment_name)


        # dividend process coeff: inflation
        _reset_params_orig()
        experiment_name = "dividend_pi"
        param_name = "p_dpi"
        param_type = :est
        param_vals = Vector{Float64}(range(-0.5, 0.5, length=N_param_vals_robustness))
        df_qe_experiment = run_experiment_robustness(
            param_name, param_type, param_vals, qe_shock, nkphopt_qe; tau_resp=tau_robustness
        )
        _df2xls_sheet!(df_qe_experiment, xls_output_fpath, "qe_exp_" * experiment_name)

        # dividend process coeff: mean reversion
        _reset_params_orig()
        experiment_name = "dividend_reversion"
        param_name = "k_d"
        param_type = :est
        param_vals = Vector{Float64}(range(0.2, 1.5, length=N_param_vals_robustness))
        df_qe_experiment = run_experiment_robustness(
            param_name, param_type, param_vals, qe_shock, nkphopt_qe; tau_resp=tau_robustness
        )
        _df2xls_sheet!(df_qe_experiment, xls_output_fpath, "qe_exp_" * experiment_name)
    end

    # save QE shock size for alternative models
    return qe_shock_size
end


#############################################################################################
# wrapper
function summarize_model_fit(model_fname::String, is_high_crisis::Bool, xls_fname::String)
    # load model and summarize fit
    # all experiments saved to xls file
    xls_fpath = OUTPUT_DIR * xls_fname * ".xlsx"
    # load model
    nkphopt = load_nkph_model_estimates(model_fname, is_high_crisis, display=true)
    # create xls file, save basic info
    df_model_summary = DataFrame(
        description = ["root_max", "qe_model_fname"],
        values = [maximum(abs.(nkphopt.nkphsolvr.root_vec)), model_fname]
    )
    XLSX.writetable(xls_fpath, overwrite=true, "summary" => df_model_summary)
    # scalar moments
    df_moment_mat = _get_scalar_moments(nkphopt)
    if size(df_moment_mat, 1) > 0
        _df2xls_sheet!(df_moment_mat, xls_fpath, "scalar_moments")
    end

    # nonscalar moments
    df_mmts_arr, title_str_arr = _get_nonscalar_moments(nkphopt)
    for (df_mmts_nonscalar, title_str) in zip(df_mmts_arr, title_str_arr)
        _df2xls_sheet!(df_mmts_nonscalar, xls_fpath, title_str)
    end

    # save alt localization and risky localization results
    sig_bl, sig_bs = estimate_sig_auc_size(nkphopt)
    tau_arr_reg = nkphopt.tau_arr_reg
    # larger set of maturities (avoid divide by zero errors)
    tau_risky_arr = Vector{Float64}(range(0.0, 30., step=0.1))
    tau_risky_arr[1] = 1e-7
    nkphsolvr = nkphopt.nkphsolvr
    J = nkphsolvr.J
    for j=5:6
        y_shock = zeros(J)
        if j==5
            y_shock[j] = sig_bs
        else 
            y_shock[j] = sig_bl
        end
        yc_resp, yc_resp_til = calc_yield_response(y_shock, tau_arr_reg, nkphsolvr)
        yc_resp_full, yc_resp_til_full = calc_yield_response(y_shock, tau_risky_arr, nkphsolvr)
        # basis points
        yc_resp = 100 .* yc_resp
        yc_resp_til = 100 .* yc_resp_til
        yc_resp_full = 100 .* yc_resp_full
        yc_resp_til_full = 100 .* yc_resp_til_full
        # alt localization
        if j==5
            alt_phreg_coeffs = nkphopt.nkphemp.df_alt_regs[:, :alt_localization_short]
        else
            alt_phreg_coeffs = nkphopt.nkphemp.df_alt_regs[:, :alt_localization_long]
        end
        # combine results in dataframe with maturities and save to excel
        alt_phreg_mat = hcat(tau_arr_reg, yc_resp, alt_phreg_coeffs)
        risky_phreg_mat = hcat(tau_risky_arr, yc_resp_full, yc_resp_til_full)

        # save
        df_alt_phreg_mat = DataFrame(alt_phreg_mat, [:tau, :b, :bhat])
        sheetname = "alt_phreg_" * string(j)
        _df2xls_sheet!(df_alt_phreg_mat, xls_fpath, sheetname)
        
        df_risky_phreg_mat = DataFrame(risky_phreg_mat, [:tau, :y, :y_til])
        sheetname = "risky_phreg_" * string(j)
        _df2xls_sheet!(df_risky_phreg_mat, xls_fpath, sheetname)
    end

    return nothing
end


#############################################################################################
#############################################################################################

println("parsing options")
parsed_args = parse_args(ARGS, s)

RUN_SUMMARY = parsed_args["run_summary"]==1
RUN_BASELINE = parsed_args["run_baseline"]==1
RUN_ROBUSTNESS = parsed_args["run_robustness"]==1
# number of times to run each robustness exercise
N_param_vals_robustness = parsed_args["N_param_vals_robustness"]



#############################################
#############################################
if RUN_SUMMARY
println("="^100); println("="^100)
println("Summarizing model fit (low crisis)")
summarize_model_fit("baseline_low_mmts", false, "low_crisis_fit")

println("="^100); println("="^100)
println("Summarizing model fit (high crisis)")
summarize_model_fit("baseline_high_mmts", true, "high_crisis_fit")
end


#############################################
#############################################
if RUN_BASELINE
println("="^100); println("="^100)
println("Running baseline experiments")

println("QE1")
# note: get size of QE1 to compute size of QT
qe1_shock_size = run_qe_experiments_all(
    "baseline_high_QE_passive", "experiments_QE1";
    run_baseline=true, run_robustness=false, run_macro=false,
    qe_shock_size=nothing, qe_shock_idx=nothing,
    mp_shock_size=-1.0, mp_shock_idx=1
)

#############################################
println("QT")
qt_shock_size = -2.0 * qe1_shock_size
run_qe_experiments_all(
    "baseline_med_QT_passive", "experiments_QT_passive_med";
    run_baseline=true, run_robustness=false, run_macro=false,
    qe_shock_size=qt_shock_size, qe_shock_idx=nothing,
    mp_shock_size=-1.0, mp_shock_idx=1
)


#############################################
# RN MP experiments
println("RN QE/MP")
run_qe_experiments_all(
    "baselineRN_high_QE_passive", "experiments_MP_QE_RN";
    run_baseline=true, run_robustness=false, run_macro=false,
    qe_shock_size=qe1_shock_size, qe_shock_idx=nothing,
    mp_shock_size=-0.5, mp_shock_idx=1
)
run_mp_experiments("baselineRN_MP_shock", "experiments_RN_MP_shock";
    mp_shock_size=-0.50, mp_shock_idx=nothing
)

end


#############################################
#############################################
if RUN_ROBUSTNESS

println("="^100); println("="^100)
println("Running QE1 robustness")
run_qe_experiments_all(
    "baseline_high_QE_passive", "experiments_QE1";
    run_baseline=false, run_robustness=true, run_macro=true,
    N_param_vals_robustness=N_param_vals_robustness, tau_robustness=10.,
    qe_shock_size=nothing, qe_shock_idx=nothing,
    mp_shock_size=-1.0, mp_shock_idx=1
)

println("="^100); println("="^100)
println("Running QE1 robustness (safe only)")
run_qe_experiments_all(
    "baseline_high_QE_passive_safeonly", "experiments_QE1_safeonly";
    run_baseline=true, run_robustness=true, run_macro=true,
    N_param_vals_robustness=N_param_vals_robustness, tau_robustness=10.,
    qe_shock_size=nothing, qe_shock_idx=nothing,
    mp_shock_size=-1.0, mp_shock_idx=1
)

println("="^100); println("="^100)
println("Running QE1 robustness (risky only)")
run_qe_experiments_all(
    "baseline_high_QE_passive_riskyonly", "experiments_QE1_riskyonly";
    run_baseline=true, run_robustness=true, run_macro=true,
    N_param_vals_robustness=N_param_vals_robustness, tau_robustness=10.,
    qe_shock_size=nothing, qe_shock_idx=nothing,
    mp_shock_size=-1.0, mp_shock_idx=1
)

end



